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BEAMS ON ELASTIC SUPPORTS AND ON 


CROSS-GIRDERS 
by 
J. M. KLITCHIEFF 


(Professor of Applied Mechanics, The Technical University, Belgrade) 


SUMMARY 


The beam on elastic supports is often considered as a beam on a continuous 
elastic medium. This method can be applied only to a large number of equally 
spaced supports having equal rigidities. The method proposed by the author, 
giving the solution in the form of rapidly converging series, is not bound by this 
restriction, and when applied to the problem of beams supported by cross-girders 
reduces it to a system of linear equations, their number being equal to the number 
of cross-girders, while the existing approximate method leads to the same number 
of simultaneous differential equations. 


Notation 


m—1 


u 
M, M’, 


2. ..., m—1 


i; 
Ri, R’, R” ;, R”, 


l 
L 
I 


J 


total number of lateral loads, or of intermediate supports, or 
of cross-girders 


distance of cross-section considered, measured from left-hand 
end of beam 


deflection of beam at cross-section considered 
bending moments in cross-section considered 


number corresponding to a particular lateral load, or inter- 
mediate support or cross-girder 


distance of cross-section denoted by i from left-hand end of 
beam 


deflection of beam at cross-section denoted by i 


lateral load, or pressure exerted by beam on support, or on 
cross-girder as support 


length of beam 
length of cross-girder 
moment of inertia of cross-section of beam 


moment of inertia of cross-section of cross-girder 
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E  Young’s modulus 
a, % spring constant of elastic support, or supports if equal 
B P/EI 
L'/EJ 
Q _ given lateral load acting on one beam 
j integers 


A, coefficients of trigonometric series representing bending 
moment in beam 


Crs C’n, ng coefficients of trigonometric series representing deflection curve 
of girder 


S, given by equation (11) 


y.8 coefficients of deflection depending on distribution of load 
and of cross-girders along beam respectively 


|n—2m| absolute value of n—2m 
6 nr/(2m) 
K,, L, see equations (20) 


1. PRELIMINARY 
The following relations are used in this paper. 


(@) If m—1 lateral loads act on a beam simply supported at the ends the bending 
moment at a distance x from the left-hand end can be expressed by the 
trigonometric series 


sin—, . (1) 
in which / denotes the span of the beam and €;, the distance of the load R, from 
the left-hand end. 


(ii) The sum 


m—1 

ime . AE nz 

in —— =sin? — cot <— (2) 
1 m 2 2m’ 


iM 
~” 


in which n and m are integers, is equal to zero for even values of n and for 
values of n which are multiples of m. For all other values of rn it is equal to 
cot nx/(2m). 
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Fig. 1. 
(iii) The sum 
1 ike ™1 f(n—k)x 


m—l inn 
sin — sin — cos 
i=1 m m i=1 m i=1 m 


in which n,k and m are integers, is equal to m for k=n+2 jm, and to —m for 
k= —n+2jm (j=1, 2, 3). It is equal to zero for all other values of k and, 
also, for values of n which are multiples of m"’. 


(iv) By differentiating three times the known series 


= 
x =cotd 

it can be proved that 


x 1 x nn ( nn 
(2jm+ny (! 1 +3 cot ~~). (4) 


2. BENDING OF A BEAM ON ELASTIC SUPPORTS 


The problem (Fig. 1) of bending of a continuous girder simply supported at 
the ends on rigid supports and having m— | intermediate elastic supports, is usually 
discussed by replacing the action of the elastic supports by the action of a continuous 
elastic medium. The use of this method is limited to the cases of many equally 
spaced supports of equal rigidities. As shown later these restrictions become 
superfluous if the expression of the deflection curve of the beam in the form of 
trigonometric series is used. 


The bending moment produced by any given concentrated or distributed load 
can be written in the form of a series 


M’= ——= A, sin —, 2 (5) 
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in which A,, are given coefficients. Denoting by u; the deflection of the beam at the 
support number i and by a; the “spring-constant” of this support, the bending 
moment produced by m—1 reactions R;=2,u; of the elastic supports becomes 
(equation (1)) 


If the deflection curve of the beam is assumed to be of the form 


u= XC, sin — (7) 
n=1 
and is introduced into equation (6), this becomes 


1 l n=1" l l 


Introducing the expressions (5), (7) and (8) in the differential equation of the 
deflection curve 
au 


=~ (M’+M") 


where El] is the rigidity of the beam, a system of linear equations can be obtained 
of the form 
2 


a 


.. nee, 
> 2, sin gin =n A,, (9) 


4 
n*cC,+ = 


k=1 


in which B=I°/EI. This system can be solved easily by successive approximations 
in every particular case. 


In the special case of equally spaced supports of equal rigidities 2, can be 
replaced by z, ¢ by li/m. If the values of relation (3) already established are 
introduced, the system of equations (9) becomes 


n*C,+ 2 BS,=n" As, 
(n— C\n—2m| — BS,=(n— 2m)? Ajn—am\, 


(n+ 2m) 2 BS, =(n+ Anam, 


(n— 4m)* C\n—am| — 2 BSy=(n— 4m)? 


il 
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in which the following notation is introduced 


Eliminating S, from the first pair of equations (10) 


n* n’ A, A \n-2m| 


and proceeding in the same way, 


Me A\n-4m 

Cin-4m| = — 


Thus, equation (11) becomes 


1 As = A\n-2m| 


Then, from the first of equations (11), using relation (4), the final equation 
determining C, can be obtained 


(a-4 [1+ (1+ )(1+3cot? 


3 A \n-2im| _$ ] (12) 
n‘ x* =1(n—2jmy j=0 +2jmy ‘ 


3. CROSS-GIRDERS SUPPORTING SEVERAL BEAMS 


The method already used can be applied in the case of several equally spaced 
and equally loaded beams of equal rigidities, supported by a longitudinal girder 
(Fig. 2). The deflection u; at the joint number i can be written in the form 


in which Q denotes the load acting on one beam, R,; the pressure exerted by the 
beam on the girder, and y and 8 are coefficients depending on the distribution of 
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it 


Fig. 2. 


the load Q along the beam and the position of the cross-girder. For example, if 
the load is distributed uniformly and the girder bisects the beams, y=5/384 and 
8=1/48. From equation (13) 


kri 
The bending moment produced by the pressures of m—1 beams is 


Introducing this into the differential equation of the deflection curve of the 
cross-girder and proceeding as in the case of elastic supports, 


oo m—1 ; m—1 


C, sin— sin — = S sin — 
BB m m B i=1 m’ 


in which B’=L*/EJ, J being the moment of inertia of the cross-section of the girder. 
Using relations (2), (3) and (4) this equation becomes 


( 2 |= cotnz/2m 
C.[ 1+ 1 +cot (1+3 cot (15) 


R; being determined by equation (14). 


Equations (15) and (14) can be applied also to the case of two symmetrically 
situated girders of equal rigidities. The coefficient 8 in this case must correspond 
to the deflection of the beam produced by two equal and symmetrically situated loads. 


4. BEAMS SUPPORTED BY SEVERAL GIRDERS 


The same method can be applied to the case of several cross-girders, and will 
be illustrated by the following numerical example. 
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Fig. 3. 


Consider (Fig. 3) the case of 11 equally spaced beams bearing uniformly 
distributed loads Q and supported by 5 equally spaced cross-girders, the moment of 
inertia of the central girder being twice as large as the moments of inertia of the 
side-girders. This might correspond, for example, to the bottom structures of a 
| large ship. 


Denote by R’; the pressures exerted by the beam number i on the girders I and 
V, by R”; the pressures exerted by it on the girders II and IV, and by 2 R”’; the 
pressure on the central girder. The beam can be considered as being acted on 
simultaneously by the given load Q, two equal reactions R’; at distances //6 from 
its ends, two equal reactions R”; at distances //3 and two equal reactions R”’; at 
distances //2. 


The deflections of the beam at its joints will be respectively 
u”,=(y" Q-B,, R’s— 8.5 . (16) 
This system of equations corresponds to equation (13). 


The numerical values of the coefficients y can be calculated from the known 
expression for the deflection of a simply supported beam produced by a uniformly 


distributed load, 


Replacing x by //6, 21/6, 31/6, respectively, we obtain 


_,_ 205 
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The values of the coefficients 8 can be calculated from a similar expression for 
the deflection produced by two equal concentrated loads at distances £>x from 
its ends 


giving 


Thus, equation (16) becomes 


14 R’,+23 R”;+26 
23 R’,+ 40 +46 


26 R’; +46 +54 


Solving this system we have 


1 118 6* , ” } 
1 


53 6* ae 


R’; 


corresponding to equation (14); or, writing 


oo 
C’, sin—— , C”, sin—, C’”’, sin —, 


L n=1 L n=1 b& 


this becomes 


1 118 6* kri 


1 100 ” kri 


” 
= 


R 
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Introducing these values in the differential equation of the deflection curve of 


the first girder, the following is obtained. 
il ; 4 
2B’ sin B 2x 118 nx 12x6 


in which 
Ss’. aC". + +...-C’ \n—2m —C’ |n—4m| — 


12 x 6* _ 2QOB' 118 nt 
C’,4 36 x* B (44C’, 46C 24 4° 


Using the notation 
26x48 x12 BN! 5g + 3 cote). 


(20) 


this can be written in the form 


. 


In the same way for the cross-girders II and III 


(22) 


—46 K,,C’,+(1+80 QB, 
18 K, L, QB. 


This system of three equations corresponds to equation (15). The solutions are 


118 + 16744 K,, + 138500 K,,? 


C.= T5155 K, +2808 K,2+616K,° 


L, QB, 
100 + 15366 K,, + 237952 K,,” 
1+155 K, +2808 K,” + 676 K,° 


_53+9048 K,, +273780 K,’ 
1+155 K, +2808 + 676 K,° 


L, QB. 
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From equations (18), (21) and (22) final expressions for the reactions R; can be 


obtained, and give the solution of the problem 


R= 1546 + 192844 K, + 79768 L, sin} 
pais... 1 +155K, +2808 Ke? +670K, "125° 

R’ Q 100 — 6* 134+ 42848 K,, + 67600 K,,” L sin 
26 1 +155 Ka + 2008 

= 2 {53 124956 K, + 35828 K,? sin } 
ants... 1+ 155K, +2808 K,?+616K, 


CONCLUSIONS 


The problem discussed contained 18 statically indeterminate reactions. 
They could be found by solving 18 equations, containing 8 unknown quantities each, 
or using the approximate method proposed by I. G. Boobnov”’. Replacing the 
pressures exerted by the beams on the cross-girders by the pressure of a continuous 
elastic medium reduces the problem to a system of 3 simultaneous differential 
equations of the fourth order. This system, in principle, presents no difficulties, 
but it can hardly be expected that the mathematical work involved can be done by 
a designer. 


The method used in this paper reduces the problem to two systems of 3 linear 
algebraic equations each. In the general case, if no symmetry exists, the number 
of equations is equal to the number of cross-girders. 
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THE ESTIMATION OF RANGE OF JET- 
PROPELLED AIRCRAFT 
by 


H. PEARSON, B.A., A.F.R.Ae.S. 
(Rolls-Royce Limited, Derby) 


INTRODUCTION 


The statement that with modern aircraft and jet propulsion engines it is essential 
to match the engine and aircraft has become a platitude, without the reasons being 
quite so well appreciated. It is well known that should the aircraft be operating 
away from its “design point,” involving throttling back of the engine, then the 
loss of efficiency with either jet or turbine airscrew engines is liable to be serious, 
and to have a large effect on the range. There is also a firm conviction that very 
high altitude is necessary in order to give a jet aircraft any range worth mentioning. 


It seems worth while to examine in a simple way on what qualities of the engine 
the range depends. The object of this simplified attack on the problem is two-fold; 
first of all it enables the factors to be readily appreciated, and second it provides rapid 
methods of comparing the qualities of two different engine types for range. 


An aircraft manufacturer or operator will always be interested in the 
performance of a particular design of aircraft fitted with particular engines, and 
this paper is not intended in any way to replace or simplify such essentially detailed 
calculations. What is intended is a simple method of estimating the maximum range 
under given conditions of any engine type applied to an aircraft of given 
characteristics. 


A large part of the paper will not be new to many interested in aircraft 
performance, but a number of the conclusions and “ready reckoner” methods, and 
the general line of argument, are certainly not generally appreciated. The methods 
are particularly applicable to jet engines but are easily applied to all types of engine. 


Notation 
Aircraft data 
D aircraft total drag (Ib. weight) 


o relative air density 


Paper received February 1950. 
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aircraft speed (m.p.h.) 

aircraft weight (Ib.) 

equivalent air speed for minimum drag 

(V/v)/o, ratio of speed to minimum drag speed at that altitude 
still-air range (miles) 

fuel capacity of aircraft (Ib. weight) 

(E + F)/W, proportion of aircraft weight allotted to power plant and fuel 
aircraft drag coefficients 

optimum lift/drag ratio 

aircraft efficiency factor 

aspect ratio 

profile drag coefficient 

wing loading (Ib. wt./ft.2)=W/S 

sea level air density=0.0765 Ib. /ft.* 

wing area (ft.?) 

profile drag coefficient of lifting surfaces 

non-lifting surface drag at unit equivalent air speed 

D3 /(SCpw), ratio of non-lifting /lifting surface drag 


engine weight (Ib.) 

engine thrust (Ib weight) 

specific fuel consumption (Ib. /hr.)/Ib. of thrust 

specific thrust extrapolated to sea lével= X /(cE) 

coefficient of variation of specific fuel consumption with forward speed 


specific fuel consumption based on power ( (Ib. /hr.)/horse-power) 


refers to conditions at the commencement of flight 
refers to minimum conditions 
refers to maximum conditions 
refers to optimum conditions 


refers to critical conditions in stratosphere 
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|. GENERAL AIRCRAFT DRAG EQUATION AND CHARACTERISTICS 


In the greater part of this paper it is assumed that the effect of Mach number 
on drag may be neglected. This is justified by the fact that it is generally accepted 
that it never pays to cruise under conditions where compressibility has increased 
the drag. Accordingly the results of this paper may be applied, subject to the 
restriction that predicted speeds are less than the critical values. The effect of this 
restriction is examined in Section 7. 


Under these conditions the aircraft drag is made up of profile and induced drag. 
D=ac V?+bW?/(c 
where D is the total drag, a and b are constants, o is the relative density, V the 
aircraft speed, and W is the aircraft weight. 


The speed of minimum drag is obtained by differentiation and equating to zero. 
acV=bW?/(c V*) 


ie. minimum drag is obtained when the profile and induced drags are the same. 
If v is put equal to the equivalent air speed for minimum drag, then 


av?=b W?/0*=4 Dain =W/QN1) 


where N, is the optimum lift/drag ratio. 
It follows that the aircraft drag may now be written in non-dimensional form 
2N,D/W=n?+1/n? . i (1) 
where n is the ratio 0 V/v or the ratio of aircraft speed to that of minimum drag 
at that altitude. 


Equation (1) is the most convenient form in which to have the basic drag 
equation. For all future purposes of this paper, an aircraft is to be specified by 
three parameters 


N, the optimum lift/drag ratio, 


W (or W,) _ the aircraft weight which will be a function of time, since fuel is being 
used up in flight. W, is therefore the weight at the start of flight, 


v (or v,) the equivalent air speed for minimum drag which again will vary 
somewhat in flight, and so v, is the value at the start of flight. 


In a propulsive system where efficiency is independent of forward speed, 
i.e. an idealised piston or airscrew engine, then since the drag times the range is 
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the work done, which for a fixed amount of fuel must be constant, it is clear that for 
a maximum range one must fly at minimum drag speed. So long as this is so the 
actual altitude of flight does not matter, so far as the aircraft is concerned. 


For the greatest endurance it is essential to fly at the minimum power, i.e. 
DV must be a minimum, 


2NiD_ = ( v 
or Ww” must be a minimum. 
This means that for maximum endurance the aircraft should be flown at the lowest 
possible altitude, and at a value of n of 1/ ¢/3. Continuous flight at this speed may 
not be possible and in any case the assumed drag equation may no longer be true so 
near the stall. 


On the other hand, for a propulsive system where specific fuel consumption 
based on thrust does not vary with speed, i.e. an idealised jet engine, the results are 
somewhat different. If the specific fuel consumption is v then the time of flight for 
a given fuel load F is F/(v D) and the range is 


R=FV/(vD) 


WJ/cR F 
Hence the maximum range is obtained for a minimum value of n+1/n’, ie. at 
n= %/3. It should be emphasised that this is the maximum range speed for a constant 
altitude when the pilot is able to adjust the throttle to control the speed. In 
general it is not advisable to operate a turbine engine at other than maximum 
continuous cruising throttle setting, in which case the altitude is no longer a free 
variable. It will be found later that for this case n= {3 is no longer the best 
cruising speed for the aircraft. 


For maximum endurance F/(v D) should be a maximum, that is, the aircraft 
flies at the minimum drag speed. 


It should be remembered that for a given value of n at which the aircraft is 
flying, the actual drag is independent of the altitude. The power required, however, 
increases with altitude. 


2. SUMMARY OF CONDITIONS GOVERNING RANGE VARIATION 
WITH ALTITUDE 


It seems worth while now to state in a compact form the reasons why range 
varies as it does for the various types of power plant. 
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Fundamentally the aircraft itself does not gain in efficiency from flying at high 
altitude. It is often stated that jet aircraft fly at high altitudes because the drag is 
low, but actually at the minimum drag speed the drag is independent of altitude. 
Thus for the aircraft itself range is uninfluenced by the altitude of flight. 


These considerations apply more or less for a piston-engined aeroplane. The 
only factors which tend to produce an optimum altitude for flight are that for a 
given aircraft there is usually some variation of airscrew efficiency with forward 
speed and a variation of the specific output of the engines with altitude. This latter 
effect is very serious above the full throttle height and generally speaking, it never 
pays to design for cruising an aircraft above the full throttle height of the engines 
because this is uneconomical in engine weight and so detracts from the 
possible fuel load. 


For jet engines the picture is very different and apt to be confusing if viewed 
from the point of view of power, but the basic reason for jet aircraft having an 
optimum range at very high altitudes is that as the forward speed increases, the 
thrust remaining constant, the power developed and therefore the efficiency, increases 
directly as the speed. A limit to the altitude is set, often by Mach number 
considerations, but also by the fact that as a constant thrust at all altitudes is 
required, bigger and bigger power plants have to be installed as the altitude is 
increased. Thus the room left for fuel is reduced and so the range eventually falls. 
The actual optimum altitude for a particular aircraft is determined by the specific 
weight of the power plant. 


The remaining sections of this paper give analytical expression to these 
arguments. The aim is to derive the range and the optimum conditions. No account 
is taken of the equally valuable cruising speed. In many cases it will pay to accept 
a reduced range at higher speed, in fact this is the selling point of jet propulsion 
power plants. However it cannot be allowed for in any simple way that has any 
real bearing on air transport, and in any case the paper gives readily usable 
expressions for this quantity. 


3. BASIC ENGINE PARAMETERS 


To propel an aircraft at any given speed and altitude a definite power is required, 
that is, either a definite size of engine is fitted, or the engine is throttled to suit these 
conditions. In this paper only maximum range conditions are considered, when 
it is undesirable to throttle. Thus it is assumed that for a given aircraft the amount 
of engine that is installed is always adjusted to meet the conditions. The only 
parameters of an engine that interest us therefore are the specific fuel consumption 
and the thrust for a given engine weight (specific thrust). 


The specific fuel consumption v for a jet engine is usually defined as the fuel 
consumption (lb. wt.) per hour per Ib. of thrust. The fuel consumption and thrust 
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being measured in the same units however, specific fuel consumption is to this 
extent independent of the units used. It is felt that some apology is necessary for 
such a definition since, unlike the specific consumption of a piston engine, it is 
by no means a measure of efficiency. It is, in fact, dimensionably an energy per 
second divided by a momentum per second, and thus has the units of velocity. The 
velocity in question is that of the aircraft. The great use of such a term is that it 
varies very little with flight conditions, and for many approximate purposes may be 
considered constant. 


The specific thrust X/E is the thrust at any specified conditions divided by the 
engine weight. It is taken as independent of engine weight E itself, although for a 
given engine type X/E may easily vary with E. Such variation within the useful 
range is normally small and is taken into consideration when the actual engine size 
required is fixed. 


Most of the quantities used are usually required in terms of the aircraft weight, 
that is E/ W is the fractional engine weight of the aircraft. For clarity it is preferable 
to leave the ratio as written rather than to give it some obscure symbol. 


In considering the relative performance of various power plants the method 
adopted is to keep the aircraft characteristics constant. This automatically means 
keeping the ratio (E+ F)/W, constant for various power plants. It is not inferred 
that this is automatically the best thing to do, for example if a piston engine were 
being replaced by a high performance jet engine it might be advantageous to increase 
the value of (E+F)/W, thus increasing the wing loading. Such a change, from the 
point of view under discussion, would constitute a new aircraft since the minimum 
drag speed would have been changed. 


4. RANGE EQUATIONS WITH ENGINE PARAMETERS 


Since (E+F)/W is to be regarded as constant (=C), the equation for range 
becomes 


This equation is reasonably satisfactory for small values of F/W=C—-—E/W, but for 
larger values, the progressive lightening of the aircraft as the fuel is consumed has 
an important effect in increasing the range. This will be taken into consideration 
later. For the moment it is easier to use the simple range formula (2). 


It is necessary to introduce into (2) the fact that it is undesirable to operate 
the engine at any other condition than maximum continuous cruising. To do this 
it is assumed that 


. ‘ > (3) 
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k is the specific thrust corresponding to sea level conditions. It assumes that the 
engine thrust is independent of forward speed and directly proportional to «. The 
former is certainly a good enough approximation, especially for the range in 
which we are interested. The second is not sufficiently good to use over a wide 
range of altitude; it is of course true in the stratosphere, and good enough elsewhere 
over a small range of altitude. Thus the value of k as given is actually the specific 
thrust extrapolated to sea level, rather than the actual sea level value. 


In the foregoing range formula there are three quantities at one’s disposal, o, 
E/W and n. For fixed throttle setting there is the second relation (3). 


Now k= = IN, (n +1] /n )} E) 
1 
oc Wnri+1/n? 
_ n E 


All parameters are now independent. 


If it is assumed that v is independent of forward speed then an optimum value 
of n= /2 is found by normal methods. It is generally accepted that for jet-engined 
aircraft the best cruising speed is 1.18 times the minimum drag speed, provided that 
the aircraft is powered suitably so that the engines are operating at the 
correct conditions 


An extremely good approximation for the variation of v with forward speed is : — 
v=Votyn, 


where y is a constant depending both on the aircraft and engine. Using this a new 
optimum value could be found for n, but the value of y is usually so small that its 
effect may be neglected. For a piston-engine power plant, the variation of v is so 
important that the estimation of range would have to be carried out graphically. 
In this case it is much better to derive the equations on the basis of specific power 
fuel consumption as is usually done for such engines. 


Using equations (4) and (5) the range for a given power plant type may now 
be plotted against altitude. The range of all types tends to increase with altitude, 
reach an optimum value, and then decrease. The reason for such an optimum is 
that the falling atmospheric pressure and constant thrust require larger and larger 
engines, which in the end swamps the benefit obtained by flying faster. 
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Attention will be concentrated first on the optimum altitudes for the different 
engine types. Differentiating equation (5) with respect to E/W and setting to zero, 
gives the optimum value of E/W 


(E/ =C/3, . ‘ (6) 


that is, at the optimum altitude the power plant weight of all jet type engines is the 
same for a given value of (E+ F) and is equal to one third of this value. In fact 
the engine weight should be one half the fuel weight. It will be seen later that 
this surprisingly simple result is somewhat modified when allowance is made for 
using up fuel; nevertheless the result is sufficiently accurate to explode a common 
belief that jet aircraft will always carry much heavier fuel loads and have 
correspondingly lighter engines than the equivalent piston-engined aircraft. This 
is only true when the jet aircraft is designed to operate at much lower altitudes 
than its optimum value. 


For example, if the design altitude were such that the air density were exactly 
double the optimum value, the value of E/W would be halved, and the fuel weight 
would be five times the engine weight. Under these conditions the range would be 
about 89 per cent. of the optimum value. Doubling the air density corresponds 
to a reduction of altitude from say, 48,000 ft. to 33,000 ft. 


The actual optimum altitude is given from equations (4) and (6), 


9 | 
Toot = 
and the optimum range is 


It will be seen that a criterion for range for a power plant is thus the value of 
Vk/v, the square root of the specific thrust extrapolated to sea level, divided by 
the specific fuel consumption under flight conditions. This seems to be a 
quick “ready reckoner” method for the power plant engineer to assess 
different engine types. 


For the aircraft the optimum range is given by the maximum value of N,? v. 
This point will be referred to later. 


5. EFFECT OF REDUCTION IN AIRCRAFT WEIGHT DUE TO 
CONSUMPTION OF FUEL 


For a long-range aircraft in which the fuel might be as high as 40 per cent. 
of the aircraft weight, the change of weight as the aircraft flies is no longer 
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negligible. It is easily allowed for by writing equation (2) in differential form 


n 
ann (=#). 


Eliminating « by means of equation (4), 


Now v is still a variable with aircraft weight, actually voc / W, so that 


+1/n?P?? w.) ( w 


and integrating between the limits of full and empty fuel load gives 


sj2,, __ _c 


This is the formula equivalent to (5) taking full advantage of progressive weight 
reduction. It assumes that the aircraft is climbing at such a rate as to keep at 
the correct altitude, and in fact the cruising speed remains constant. 


Again an optimum value of the range is obtained for cruising at n= //2. 
Fig. 1 shows plots of equation (9) made in such a manner as to be universal for all 
aircraft and power plants. This is done by plotting “non-dimensional” range 
v, Vk/¥) against a sort of “non-dimensional” altitude 1/(¢k N,). Thus 
it is only necessary to substitute actual values of N; and wv, for a particular aircraft, 
and /k/v for a power plant to obtain the range plot for an actual combination. 
It will be seen that the range rises to a maximum value and then falls again, as 
explained earlier. 


By differentiating (9) with respect to E/W and setting to zero, the condition 
for optimum range is :— 


The solution to this for various values of C=(E+F)/W, is plotted on Fig. 2, 
compared with the simple solution E/W,=C/3 derived earlier. It will be seen 
that an optimum value of E/W, is reached after which it declines rapidly. It is 
scarcely worth while at this stage to investigate the behaviour of an aircraft of 
which less than 15 per cent. consists of structure, payload and so forth. 
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o kN, 
Fig. 1. 
Variation of range with altitude. 


It will be seen that over a considerable range a fair value of E/W, to maintain 
optimum conditions is 


E/W,=0.284 C, 
in which case the approximate range is given by 


/k 
éi 
Rm =0.33Ni, 


which is the counterpart of equation (8). Again the determining factor from the 
power plant point of view is /k/v and for the aircraft the value of N,?v,. The 
very fast rise in range with increasing values of C may also be appreciated 
from Fig. 3. 
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Fig. 2. 


Variation of power plant weight with (engine plus fuel) weight for maximum range conditions. 


6 OPTIMUM DESIGN OF AIRCRAFT FOR JET PROPULSION 


It has been agreed previously that the only noticeable advantage in flying 
high is that the increased economical cruising speed is more suitable for jet engines. 
The beneficial effects of the decreased temperature on the engine performance are 
of a secondary nature. It has been suggested that if an aircraft were designed to 
have a very high minimum drag speed at low altitude a similar result would be 
obtained. A high minimum drag equivalent air speed could only be obtained by 
reducing the wing area and thus increasing markedly the wing loading. Apart from 
the feasibility of this, however, the previous range formule (8) and (11) show that 
the only criteria of range are the specific thrust, the specific fuel consumption, both 
of which are not in the least improved by such a step, and the value allotted to 
(E+F)/W, as the factors affected by the engine, and for the aircraft the 
value of N,? v. 


Thus the problem of optimum aircraft design for range with jet engines is to 
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Fig. 3. 
Variation of maximum range and optimum altitude with (E + F)/ Wo. 


make the value of N,? v a maximum. These quantities are given in terms of aircraft 


data by 


Ww 
7 2 


For the present purpose the aspect ratio is regarded as a constant and variations of 
wing area are considered. Now the aircraft body will presumably have to remain 
a fixed size as the wings are reduced, indeed the body would probably have to be 
made bigger to contain the fuel that would have been contained in the wing, but 
the body will be regarded as of fixed size and therefore as. of fixed drag D, (assessed 
at fixed speed). 


Hence Coo = Cow + D;/S > 
where Cpw is the profile drag coefficient of the wing. 
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Thus the variable term in the expression for N;? v becomes 


W and Cp being regarded as constant, this may be differentiated with respect to 
S and an optimum value is found when 


1 Ds 
S Cow 

ie. the body drag coefficient (expressed on the wing area) should be the same as 
the wing drag coefficient. Naturally this does not take account of such things as 
interference drag, but probably this ought to be included as body drag. Obviously 
the wing drag means the drag of all lifting surfaces. 


The condition required seems to be approximately met in most modern aircraft, 
as the maximum of the expression / x/(1+ x) is not very sensitive to values of x. 
For example at the optimum value of x=1 it is 0.5, whereas at x=2 and } it is 
still greater than 0.47. On the other hand, to make an aircraft have a very high 
minimum drag speed at low altitude requires a very high value of x=Dz/(S Cow), 
approximately x=29 to make the minimum drag equivalent air speed 500 m.p.h. 
Although the jet engine would be operating efficiently, the range would have been 
reduced to 36 per cent. of the optimum value because of the high body drag. 


The fact that jet propulsion exercises an optimum on the value of N,? v is 
novel and does not occur for a constant power specific engine. The equivalent 
formula to (5) for such an engine is 


_2N. 1 


where v, is the power specific fuel consumption. The conditions for an optimum 
range therefore are, cruising at minimum drag speed n=1, at the minimum value 
of E/W and therefore at sea level, assuming the specific power does not change 
much below full throttle height, and for the aircraft a maximum value of N1, i.e. 
as low a profile drag coefficient as possible. 


These conclusions are vitiated in practice, because of the variation of airscrew 
efficiency, and so forth. 


7. EFFECT OF LIMITING MACH NUMBER IN THE STRATOSPHERE 


It will be found that for most modern aircraft powered with jet engines, the 
optimum altitude for range is well above the tropopause and also, the optimum 
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speed corresponding to this altitude is above the critical Mach number of the aircraft. 
This places a limit to the possible range below the optimum figure discussed 
previously. If it is accepted therefore that there is now a constant critical aircraft 
velocity V.. at which the aircraft must be flown, the conditions for optimum range 
become somewhat different, and much simpler. 


On a simple basis 


R=(F/v\(V./D) 

so that 


Now V., n and E/W are not independent. In fact 


1_,E 2N; 

Wr'+1/n? 
and 

nv 

_2N,_ (1+ 1/n*) By 
Hence x= ~ ] 


(13) 


Vel 2NiC 


There is no obvious general optimum value of mn. It may be found in any 
specific case and it is clear that in the case of the critical Mach number limiting 
the altitude for greatest range, the greatest range will be obtained by flying at the 
critical Mach number at such an altitude that n lies between 1.18 and 1.00, the 
actual value being determined by the optimum of equation (13). 


In practice the assumption of a constant critical Mach number is an over- 
simplification, since it is likely to be a fast-changing function of the lift coeficient of 
the wings. As altitude is increased at constant speed, the lift coefficient increases 
and the limiting Mach number will fall. Such considerations are easy to take into 
consideration when the exact variation is known. 


8. PHYSICAL SIGNIFICANCE OF PREVIOUS RESULTS 


It seems worth while to discuss the results summarised in Fig. 1 in terms of 
real values applied to real aircraft. If a value of N:, the optimum lift/drag ratio, 
is taken equal to 18 and an indicated air speed for minimum drag of 180 m.p.h. 
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typical conditions will be represented. For modern jet power plants figures of k=3 
are not unreasonable. For an aircraft made up of 50 per cent. engines and fuel, 
Fig. 1 then shows that the optimum range will be attained at s=0.137 or an altitude 
of above 52,000 ft. at the start of the cruise. Even if (E+F)/W is 40 per cent. 
then this altitude is still 48,000 ft. and at the end of the flight when the fuel is 
finished it would have risen to 55,000 ft. 


It is evident therefore that to obtain the maximum range from modern jet 
engines, aircraft will have to fly very high, higher perhaps than may be permissible 
for civil use. However, as pointed out in an earlier section, the loss in flying lower 
than the optimum altitude is not too serious. The shape of the curves in Fig. 1 
shows that o may be increased greatly without a large fall-off in range. Thus, if 
instead of starting at 48,000 ft. in the foregoing example 33,000 ft. is chosen (where 
a is exactly double), the range becomes 89 per cent. of the maximum value. At 
25,000 ft. the range would be only 80 per cent. of the maximum. 


The actual maximum still-air ranges for the cases considered work out at 
3,320 miles for the 50 per cent. engine+fuel case and 2,260 miles for the 
40 per cent. case. 


The fact that the optimum altitude is so high has an important bearing on the 
increase in range for very light engines that may be developed in the future. Thus 
if one can fly at the optimum altitude, then the range has been shown to be 
proportional to /k/¥, i.e. if the value of k were increased from 3 to say, 10, then 
the range would be increased by a factor 1.81. If a limit is set to say, 
50,000 ft. as the initial cruising altitude, then this factor will be reduced to 1.3. 
Thus the returns in increased range from making very light engines will be 
considerably reduced by limitation of altitude. 


This limitation in altitude will be considerable, even for military aircraft, for 
apart from difficulties of pressurisation, and so on, the aircraft will be restricted 
by its limiting Mach number. Thus in the example quoted, supposing the limiting 
Mach number were as high as 0.85, flight would be restricted to 560 m.p.h. which 
would be attained under normal conditions of cruising, at a little above 50,000 ft. 
It would then pay to go somewhat higher at the same speed, but the returns would 
not be great. 


Thus from the point of view of increasing the range of more or less conventional 
aircraft, it would seem that it would be more profitable to concentrate improvements 
in engines on reducing the specific fuel consumption, while retaining present 
thrust/weight ratios, rather than to go for much lighter engines. 
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9. OTHER FACTORS 


In this paper conditions for maximum range only have been calculated. It is 
clear that in practice other considerations may govern the choice of power plant 
loading, for example take-off or emergency climb. The aircraft may then be 
effectively over-engined and the engine would have to be operated throttled back. 
Such considerations come outside the scope of such a general treatment as this and 
belong to the detailed calculation of performance of a specific aircraft, which it is 
not the purpose of this paper to cover. It is sufficient to say here that in general, 
for long-range jet aircraft powered to attain approximately their optimum altitude, 
there is a large amount of take-off power, and this problem should not therefore 
be met, except at high altitude aerodromes in hot climates. 


Also, in this paper no allowances have been made for fuel used in take-off, 
climb, descent, stand-off and so on. Here again these factors belong to the specific 
aircraft calculations. Their qualitative effect is not difficult to see, however; take-off 
and stand-off as well as reserve fuel represent a constant weight of fuel which cannot 
be allowed as fuel for range, and should be subtracted from the fuel load before 
the calculations are made. Climb and descent tend to cancel one another out, but 
not completely. For long range they do not play a large part in affecting the optimum 
height for flight, but on shorter range aircraft (low (E+ F)/W) they may lower the 
optimum flying height considerably. 
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ON NON-STEADY MOTION OF SLENDER BODIES 
S by 


| JOHN W. MILES 
(Associate Professor of Engineering, University of California, Los Angeles) 


SUMMARY 


Following the original work of Munk, Jones, and Ward for steady flow, a 
solution is given for unsteady, supersonic flow over very slender bodies of revolution 
and wings. The results are subject to the restriction (M?— 1)6* log. [(M?— 1) 6] 
i << where 4 is the slenderness ratio, M is the Mach number and, in addition, 
to the usual restrictions imposed by linearisation. As examples, the lifts and pitching 
moments on flat wings and bodies of revolution executing pitching motion and the 
damping moment on a rolling wing are calculated. 


It is shown that the order of approximation is consistent with the limitations 
already imposed by linearisation, at least in supersonic flow, where no Kutta 
condition is required. 


Notation 
A wing aspect ratio (see equation (28) ) 
C, lift coefficient 
Cy pitching moment coefficient 
C, rolling moment coefficient 
D /Dt differential operator with respect to time 
I, see equation (24) 
L siift (transverse) force 
Mach number 
R _ radius of cross section of body of revolution 
S area of wing plan form 
S, base area for body of revolution 
S(x) area of cross section of body of revolution 


U _ free stream velocity 
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a fraction of body length of axis of pitch aft of nose 
b span at x 

c sonic velocity in free stream 

k reduced frequency (see equation (9)) 

l body length 
m apparent mass of cross section 

p angular velocity in roll 

r polar co-ordinate 

t time 

w  downwash, i.e. velocity transverse to free stream 

x, y,z Cartesian co-ordinates 

a amplitude of pitch displacement (see equation (19) ) 
B (M?-1)} 

y  semi-apex angle of triangular wing 

8 slenderness ratio (maximum width/length) 

6 polar co-ordinate 

p density 

@ velocity potential 

rt dimensionless time (Ut/I 

® angular frequency 


1. INTRODUCTION 


The solutions for steady flow past slender bodies of revolution and narrow 
wings were first given by Munk") and Jones’) respectively. More recently, a 
critical treatment has been given by Ward"). An extension of their work to the 
problem of unsteady flow is presented here, treating only the lifting forces. While 
the results of the slender body analysis are independent, in a restricted sense, of 
Mach number, the supersonic case will be treated primarily since the effect of the 
wake may then be disregarded. 


It is assumed that the body lies in the immediate vicinity of the x axis, its 
most forward point projecting on the origin and its length being /. The free 
stream of velocity U is directed along the positive x axis, and the theory is 
linearised in the usual manner. The velocity potential then satisfies the linearised 
equation 


V*o(x, y, y, z, t) (1) 
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where c is the sonic velocity in the undisturbed medium. It is convenient to assume 
the co-ordinates (x, y, z) to be dimensionless, with / as the unit of length, and to 
introduce the dimensionless time 


Expanding the Laplacian and introducing 7, equation (1) reads 
Dyy + = + +M9--, . (3) 
where M=U/c, B=(M?-l)y. . F (4) 


The slender body approximation in steady flow assumes the flow to be 
essentially two-dimensional at every cross section, thereby implicitly neglecting 
the term 879,,. In the case of a steady flow past a body of revolution, Lighthill'*’ 
has shown that this is the best approximation which is consistently allowed within 
the framework of the linearised theory. This approximation naturally breaks down 
for sufficiently large values of 7, i.e. in the hypersonic régime but, since the 
linearised theory itself is invalid for large M, this does not represent an entirely 
new limitation. More precisely, if 5 is the slenderness ratio, the slender body 
approximation (in steady flow) requires 


whereas the hypersonic limitation requires Mé to be small, so that the two 
restrictions are not exactly equivalent. For a more critical discussion of the slender 
body approximation in steady flow, the reader is referred to the work of 
Lighthill ©) and Ward“), 


Turning to the case of unsteady flow, reference is made to a general discussion 
of the linearising process for the two-dimensional case which has been given recently 
by Lin, Reissner, and Tsien’). They conclude that necessary conditions for the 
validity of the linearised theory are 


where 6 is the thickness or amplitude of oscillation, whichever is larger, referred 
to 1, and 1/k is a characteristic time, cf. equation (9). They also conclude that 
if M is not large and either |1—M| or k >> & the linearised theory for unsteady 
flow is valid even in the transonic range. In this connection £8 is the important 
parameter in the slender body approximation in steady flow, cf. equation (5), 
whereas Ward"*) evidently assumes 8=O/(1) and uses 5 alone. 
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No attempt will be made to discuss here the sufficiency of the foregoing 
conditions for linearisation of the three-dimensional, unsteady flow problem, other 
than to point out that, in the case of the body of revolution, there is the additional 
complication of possible singularities on the axis, as a result of which logarithmic 
terms are introduced in the order of magnitude approximations. However, a very 
cursory examination of the analysis of Ref. 7 makes it evident that the restrictions 
(6) and (7) also will be necessary in the three-dimensional case. Accordingly a 
two-dimensional cross flow will be assumed and the result examined a posteriori, 
in the light of the restrictions. 


2. THE SLENDER BODY SOLUTION 


The transverse oscillations, exhibiting the harmonic time dependence exp (iw), 
of a pointed body of revolution, described in cylindrical polar co-ordinates by 
r=R (x), will be considered. The (linearised) boundary condition, to be satisfied at 
the mean position, is given by 


r(x, r, 9, T)| = —w(x)cos@exp(ikr) ® 
k=ol/U, . ‘ (9) 


where k is the reduced frequency parameter, and w(x) is the prescribed downwash. 
There is also the implicit condition that the disturbance vanish at infinity in the 
proper manner (the Sommerfeld “radiation” condition). 


If two-dimensional flow is now assumed at any cross section by neglecting the 
x derivatives in equation (3), 
+ (10) 
Drr r = rr. . . 
A solution to equations (8) and (10), which represents an outgoing disturbance at 
infinity, is given by 


w (x) H,.2) (KMr) cos 6 exp (ikr) 
KM [Hi (KMr) 


(x, r, 6, r)=—- (11) 


where H,,.) is Hankel’s function of the second kind and order one. In the vicinity 
of the body, where r~R=O (8), equation (11) reduces to 


(x, r, 9, rT) =w (x) R? (x) exp (ikr) { 1+O[(KM8) log. (kM8)}} (12) 
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Moreover, 
[ 2M*o 
O ] (M83) ‘ (13) 


It follows that, in view of the restrictions of equations (6), (7), the approximate 
solution of equation (12) is all that is justified by the linearised theory. 


While the foregoing discussion has been limited to a body of revolution, a more 
zeneral analysis could have been given, following that of Ward. Thus it is found 
that the linearised flow about any slender body may be approximated as a solution 
to Laplace’s equation provided that, in addition to the restrictions imposed 
by linearisation, both the maximum angle between any tangent plane to the body 
and the free stream direction and the rate of change of this angle are O (8) and the 
curvature of any section is O(1/d), where d is the maximum section diameter, at 
all points where the section is convex outwards (there being no similar restriction at 
points where the section is concave outwards). The last restriction guarantees 
small disturbance velocities, relative to U. If it is not satisfied, as in the case of 
thin flat wings, it can no longer be asserted that equation (12) correctly describes 
the order of the approximation, but this difficulty is present whenever the linearised 
theory is applied to thin wings, and the resultant errors in integrated results (such 
as lifts and moments) generally are acceptable. 


In considering bodies of arbitrary cross section, generally it is necessary to 
consider the interaction of the axial flow, which will introduce a logarithmic term 
in the potential, and the cross flow, as done by Ward’. In the special case of a 
body of revolution the axial flow has no net effect on the transverse forces, by 
virtue of symmetry, although naturally it must be included in the calculation of the 
drag. Thus, the concept of virtual mass may be introduced in the calculation of 
the transverse forces acting on a body of revolution and, also, on a slender, flat 
wing, since the axial flow introduces no disturbance in the latter case. (The 
application of the virtual mass concept to such problems is due to Jones’’.) In the 
case of a winged body of revolution, on the other hand, the axial flow due to wing- 
body interference must be considered, and the concept of virtual mass is no longer 
sufficient for a complete analysis. 


Since the solution to equation (3) for unsteady flow is identical with that for 
steady flow, the difference between the two cases manifests itself entirely in the 
linearised Bernouilli equation, namely :— 


pressure = p> — pU (,+ike). . 
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Thus, if m (x) is the virtual mass of a cross section at x, the transverse (lift) force 
may be calculated directly from Newton’s law, namely : — 


2 [m (x)w (x)exp (ikr)] 


=u(2 +-ik) [mm (x) w Jexp (kr) (16) 


3. BODY OF REVOLUTION 


The potential for this problem has already been obtained, cf. equation (12), 
and equation (16) is applied directly. The virtual mass of a circular cross section 
is given by 


m (x)= pS (x)= 7pR? (x) . : 
whence equation (16) reads 


w (x) 
U 


@)=pur(2 +ik)[ se Jexp Gk). 


The case of principal interest is that of a pitching oscillation of amplitude 
a about x=a, for which 


w(x)=Ua[l+ik(x-a)] . ‘ (19) 


If the lift and moment (with respect to x=a) coefficients are defiined by 
1 
0 


1 


Gkr)]-* | (@—2) L’ (21) 


0 


where S, = S(1) is the base area, equations (18) and (19) yield 


Cue=2 [—(1 —a—I,)—ik +k? . (23) 
1 
L= “5, & = é (24) 


0 
It is of particular interest that the imaginary part of Cy., i.e. the damping moment, 
vanishes for an axis through the base (a= 1). 
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4. LOW ASPECT RATIO WINGS 


In this section low aspect ratio, pointed wings, are considered, of span b (x), 
subject to the restriction 


b’ (x)>0. ‘ . (25) 


In the more general case, where equation (25) is not invoked, it would be 
necessary to impose the boundary condition 


¢ (x, y, O+)=exp [ik (x, y,04]) (26) 


at points in the wake aft of the maximum cross section, x, (y) denoting the trailing 
edge. Equation (26) ensures continuity of the pressure across the wake, even 
though the potential itself is discontinuous there. Effectively the imposition of 
equation (26) spoils the assumption of a two-dimensional cross flow, since the 
boundary conditions at sections aft of the maximum span are no longer independent. 


In the case of steady flow over a flat wing with constant incidence, Jones‘? 
and Ward) have shown that this boundary condition is satisfied by making the 
cross flow about all cross sections aft of the maximum span identical, so that no 
lift force arises (i.e. D/Dt=Ud/dx=0). Unfortunately this neat state of affairs 
does not prevail in the case under consideration, both due to the phase factor in 
equation (26) and the dependence of the local incidence, w (x)/U, on x. 


Despite the restriction of equation (25), continuity of the pressure also would 
have to be required at the trailing edge if the flow were subsonic, in satisfaction of 
the so-called Kutta condition. From the mathematical viewpoint, the Kutta 
condition, as applied in conventional aerofoil theory, is essentially a uniqueness 
condition and is required because of the existence of eigensolutions to the boundary 
value problem. In the case of the slender body solution, however, the solution for 
the two-dimensional cross flow is already unique for the boundary conditions stated, 
and no eigensolutions exist. The answer to the apparent paradox lies in the 
violation of the initial assumptions on which the slender body approximation rests, 
for if the potential is discontinuous at the trailing edge its x derivatives most 
certainly cannot be assumed small. Thus, while it has sometimes been stated that 
the slender body approximation is valid for subsonic, as well as supersonic, flow, 
this may be asserted with certainty only in those special cases* where the Kutta 
condition is satisfied automatically. While it may furnish acceptable results in 
other cases, the order of the error can be established only by a more thorough 
investigation or by comparison with experiment. Similar objections may be 


*Such as the elliptic wing used by Jones as a basis of comparison between the slender body 
approximation and a known exact result. 
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levelled at the application of the slender body approximation to the analysis of 
subsonic flow about a body of revolution having a base of finite area, but in this 
case no satisfactory theory exists to describe the actual phenomena, which must be 
heavily influenced by viscous effects. Accordingly this investigation is restricted 
to supersonic flow over wings with monotonically increasing span (which implies a 
straight trailing edge, transverse to the flow), so that no section is influenced by 
events occurring downstream. 


Returning to the special category of wings described by equation (25), the 
approach of equation (16) is used. The virtual mass of a flat plate is identical with 
that of the circumscribed, circular cylinder, whence 


If the lift and moment coefficients are referred to the plan form area S, and the 
aspect ratio is introduced 


where b, is the maximum span (in this case, the span at the trailing edge), the 
results may be taken over from those for the body of revolution by multiplying 
the latter by («A /4) and replacing S (x) by b? (x) in equation (24). 


The simplest, pointed, low aspect ratio wing, and perhaps the only one of any 
practical importance, is the triangular wing, for which 


In this case 
1 +ik -a) +k? (3 


5. NON-PLANAR MOTION 


Consider now the more general case, where the motion of the wing is not 
plane. The problem is one in potential flow (because of the slender body 
approximation) and, in principle, is straightforward. 
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y+iz= rele +h)? ele 
Fig. 1. 


Mapping of wing section on circle. 


A section of the wing is shown in Fig. 1. The (linearised) boundary 
conditions are 


y, Z)| <0. = —wx, |y|<4b(x) 


o(x, y, 0)=0, |y|>4b@) A (3) 


y,z)—> 0, |y z|—> oo . (34) 


while ¢ is required to satisfy Laplace’s equation 
Pyyt%2=0. . : . (35) 
The harmonic time dependence factor has been suppressed. 


The desired solution is obtained by mapping the wing (z=0, | y| << b/2) on 
the circle r=(b/4) with the aid of the Joukowsky transformation 


2 


r > 


as shown in Fig. 1. 


f z 
| 
r 
@=0 
| (x) +/>(x) y 
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Proceeding in the usual manner, a solution to equations (32) to (35) is given by 


b | ( b 7) si 
9(x, 20084, z w X, COS sinngsingdg . (37a) 


b 1+(4r/ —2 (4r/b) cos (6+¢) 
~ On Jw( x. $008 sin ¢ log. 
(37b) 


On the aerofoil (r= 5/4), equation (37) reduces to 


(x, 6, 0+) = >» w (x, cos ) sin ng sin ¢ d2 (38a) 
1 


< 
0 
© 


x, 5 COS ¢ sin log, dz. . (38b) 


For the lift and rolling moment on the wing section at x, the following 
integrals are required 


(2x) 
(x)= o(x,y,0+)dy . . (39a) 
—4b(z) 
= PO)! 0+) sind dd : . (39b) 
0 
0 
4b (x) 
(x)= | o(x, y, 0+) ydy. , . (40a) 
— 4b (z) 
=7 x, 6, 0-+ ) sin cos . (40b) 
0 
| sin 26 sin 6 dé . (40c) 
0 


As an example the damping due to roll (angular velocity p) is considered. 
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The downwash is given by 
w (x, y)= py =p [b (x)/2] cos . (41) 


and the rolling moment coefficient by 


1 (zx) 

C.= [4pU?Sb,]- | ae [-2p(a. 0+)]ydy . (42a) 
6 (2) 
0 


Substituting the distribution (41) in (40c) and the latter in equation (425), 


Cy» = sit | . (43) 
2U 
The real and (negative) imaginary parts correspond to the damping and inertia 
respectively. For a triangular wing of semi-apex angle y, the aspect ratio A=4tany, 
and equation (43) then gives —(=/8)tany for the damping moment coefficient. The 


exact (linearised) result for the damping in roll of a triangular wing in supersonic 
flow has been calculated by Robinson“*) and, in the limit of small y, yields 


C,= tan . (44) 


in agreement with the result above for small y. The comparison reveals that the 
accuracy of the slender body approximation in predicting integrated force coefficients, 
relative to the three-dimensional, linearised theory, is correctly indicated by equation 
(5), despite the violation of the restriction on curvature of convex sections. Similar 
comparisons also have been made for the results (30), (31) with the limiting forms 
of the complete longitudinal, stability derivatives®’ (Robinson) gives only the 
quasi-stationary results, neglecting the time derivatives in the Bernoulli equation and 
the continuity equation). 
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SUMMARY 


The paper is concerned with the two-dimensional, steady, irrotational, 
isentropic flow of a perfect gas in a wind tunnel nozzle which is found to produce 
a flow in the test rhombus deviating slightly from the desired uniform flow. 


The minimum corrections are derived that must be applied to the liners in 
order to produce a uniform flow in the test rhombus. If the uncorrected nozzle 
produces a flow of uniform direction, measurement of the pressure on the axis, 
in the test rhombus, suffices to determine these corrections (Section 5). If not, 
further pressure measurements are required (Section 6). A simple test is indicated 
for checking whether the flow stream direction is uniform (Section 6). 


The method cannot be used to correct for deviations from a_ two- 
dimensional flow. 


Notation 


x,y Cartesian co-ordinates in the flow plane, with x measured in the 
downstream direction of the axis 


pressure 
density 
velocity magnitude 

Mach angle 

a, critical speed of sound 


t | 
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stream direction 
a2 
6-t 
h., hs length parameters of Mach line pattern (see Section 3) 
H. 
Hs (Agsin2uy 
H;° _ value of H; at G (Fig. 1) 
M.,M_ plus and minus Mach lines (Fig. 2) 
S streamline (Fig. 2) 
8 correction of local inclination of liner 
N 
r _ length measured on straight minus Mach line from its point of inter- 
section with the plus Mach line BA (Fig. 1) 
r, value of r on AC (Fig. 1) 
xt 
f(#) { [(y—cos 2u) (sin secu } 
y ratio of specific heats 


The foregoing symbols refer to the theoretical field of flow on which the design 
of the nozzle is based. 


™ refers to the measured field of flow in the actual (uncorrected) nozzle. 
’ denotes differences between measured and theoretical quantities. 
* refers to the field of flow in the corrected nozzle. 


” denotes differences between quantities referring to the corrected nozzle and 
those referring to the theoretical field of flow. 


1. INTRODUCTION 


Several effects in wind tunnel nozzles for two-dimensional, supersonic flow 
cannot be predicted theoretically, at present, with great accuracy. Sometimes, 
therefore, the flow obtained in the test rhombus is not as uniform as intended, and it 
may be desired to make use of the measured distribution of pressure in the working 
section in order to determine the slight corrections to the liners needed to reduce 
the deviations from uniformity of the flow below the level of experimental 
uncertainty. It is not usually necessary that the mean Mach number in the test 
rhombus should be exactly equal to that aimed at in the original design. 


It is assumed that a small correction of the curvature of the downstream part 
of the liner will not alter the displacement thickness of the boundary layer by a 
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Fig. 1. 


significant amount. It is also assumed that the flow in the uncorrected nozzle is 
indeed two-dimensional. The problem may then be treated by the theory of the 
two-dimensional, isentropic, irrotational, steady flow of a perfect gas. 


A solution of the problem on Linear Theory has been given by Ward"), a 
solution on Exact Theory by Harrop and Bright’); the latter is not complete, and 
this may have been the reason why these authors have found it necessary to apply 
their method of redesign repeatedly in order to obtain the desired uniformity of the 
flow in practice. Carriére'*’ also claims to have obtained a solution. 


The problem belongs to a class which has been treated, mainly by Guderley“, 
under the name of “ Neighbouring Solutions.” Two fields of flow are compared in 
which the respective velocity components, at corresponding points, differ by 
amounts small enough for the squares and products of such differences to be 
negligible. In the following a notation developed by one of the authors) is 
employed that yields results more simply and more directly, in the case of two- 
dimensional flow, than Guderley’s notation. 

The correction is expressed in terms of the values on the axis of the 
characteristic co-ordinates of the field of flow in the uncorrected nozzle (equations 
(*3), (14), (20), and (21), (22), (23), which may be found from a minimum of 
pressure measurements. A procedure for the evaluation of the correction is 
indicated in Sections 5 and 6. 


2. MATHEMATICAL FORMULATION 


Let A be the apex of the test rhombus on the axis of the nozzle (Fig. 1; the 
direction of the flow is from left to right), and let AB and AC be the Mach lines 
through A, and BC the segment of the liner between them. In a perfect nozzle the 
region ABC is occupied by a simple wave, since it is adjacent to a region of uniform 
flow, i.e. the test rhombus, and the Mach lines belonging to the same family as AC 
(“minus ” Mach lines) are straight in this region. 


197 


L 
B 
SAV 
4 
G R 
¢ 
U 
A Q 
| 

| 
n 

= 


R. E. MEYER AND M. HOLT 


The least part of the liner of an imperfect nozzle that must be corrected in 
order to produce a uniform flow throughout the test rhombus is the segment BC 
(and the segment placed symmetrically to it with respect to the axis), by virtue of 
the Uniqueness Theorem“. If correction is confined to it, the field upstream of, 
and including, the Mach line BA will not be affected by the correction, and therefore 
a uniform flow must be aimed at, in the test rhombus of the corrected nozzle, with 
that Mach number and that stream direction which obtain at A in the uncorrected 
nozzle.f 


For the calculation of the nozzle correction three fields of flow must be 
distinguished, differing only little from one another, namely (i) the theoretical 
solution on the basis of which the nozzle was first designed, (ii) the field of flow 
obtained in the uncorrected nozzle, and (iii) the field of flow which it is desired to 
produce in the corrected nozzle. Quantities referring to the second (measured) 
field will be distinguished in the following by an index “m” from those referring 
to the first field, and quantities referring to the third field will be distinguished by 
an asterisk. The difference between the value, at any point, of quantities in the 
second and the first field will be denoted by a dash, and the difference between 
values in the third and the first field by two dashes (for example, 6"=6+¢, 
6*=6+06”). It will be assumed that squares and products of dashed quantities 
may be neglected. 


If x and y denote Cartesian co-ordinates in the flow plane, with x measured 
in the downstream direction of the axis, and 6 denotes the angle the stream direction 
makes with the direction of x increasing, the slope of the desired liner is 


dy* /dx=tan 6* = tan (6" +3), 
with -@’, . : (1) 
and to the first order the liner correction is 
x 
*p 


where the integral has to be evaluated along the liner. 


TIf it is desired to correct so as to produce a uniform flow with the exact Mach number and 
stream direction of the original nozzle design, it becomes necessary to extend the correction 
to a part of the liner upstream of B. To compute such a correction the characteristic equations 
(18) and (19) of the perturbation problem have to be integrated step by step, by a method 
analogous to Massau’s method®, in a region upstream of BA. 


198 


as 


3. 

= 

ST 

Ww 

I 

t 


on 
ns 


CORRECTION OF WIND TUNNEL NOZZLES 


3. EQUATIONS FOR THE CORRECTION TERMS 


The differential equations governing the variation of 6 and 6” may be derived 
as follows. 


where g denotes the velocity magnitude, » the Mach angle, and a, the critical 
speed of sound. By virtue of Bernoulli’s equation, 


qdq +(1/p)dp=0, 


where p is the density and p the pressure, and we have 


If a=6+4f, 
(5) 

B=6-, 


then ‘® 


dy/dx=tan(@-) and z=constant on any “plus” Mach line, 
(6) 


dy/dx=tan (6+) and 8=constant on any “ minus” Mach line. 


As in Ref. 5 the positive direction of a Mach line is defined as that making an 
acute angle with the stream direction and the elements of length are denoted in 
the positive directions of the minus and plus Mach lines by h.dz and hsdB 
respectively. Let 


Now compare the Mach line patterns of any two fields of flow differing only 
by small quantities (the notation already introduced for the first and second field 
may be used). On the minus Mach lines of the first field d8 =0-by (5), and (Fig. 2) 


ha"dB&™ sin sin 2u = h.da (0 
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to the first order. Since the difference between As™ and hy is smallt, it follows 
that, to the first order, 


(using (6)). Similarly, on the plus Mach lines (Fig. 2), 
h.™da™ sin 2u= —hsdB (6’ — nv’), 
and hence, 
da’ —hsH. (6 —p’). 
Only one derivative occurs in each of the equations (8) and (9) and they are therefore 
characteristic equations. Hence, the Mach lines of the first field are also the 
characteristics in the physical plane of the solution for the dashed quantitiest. 
In a region of uniform flow H, and Hz both vanish identically’) and hence, 


In a simple wave region where the minus Mach lines are straight H. vanishes 
identically) and hence, 


The variation of quantities with two dashes is governed by analogous equations. 


+This may be seen as follows. By (5), 


m 
=hg™ cos(O™—w™)+ ... 


and the difference between 6™—y»™ and 6—u is also small. 


{This fact was first noticed by Guderley™ and a direct proof has been indicated by Carriére®). 
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Fig. 2. 


4. FORMULA AND PROCEDURE 


If L is any point on the liner, between B and C (Fig. 1), and G the point on 
the Mach line BA lying on the same minus Mach line as L, then (2) may be written 


*B 


since @ is constant on the straight Mach line GL, and by (1) and (5), 
This may be expressed in terms of the values of 2’ on the axis as follows. 


By (10) and (11), 2’ and 2” are constant along the plus Mach line through L 
and hence, 2’,=2’g and a”,=2”9, where Q is the point on the axis lying on the 
plus Mach line through L (Fig. 1). Moreover, «” is constant on the axis and equal 
to the value of 2’ at A, since the correction does not affect the velocity components 
at that point, so 
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Let U, U, and V, V, be respecively the points of intersection of the Mach line 
BA and of an arbitrary other plus Mach line in the simple-wave region ABC with 
any two straight minus Mach lines at a small distance apart (Fig. 3). Let O be 
the common point of these minus Mach lines, and let UU, and VV, be 
perpendicular to UV. Then, if 8=8, on UV and B=8,+d8 on U, V,, 


VV,=VV, sin 2u=hgd8 sin 2u—VO d (0 + 1), 


0 


Fig. 3. 
to the first order, and by (5), (3), and Bernoulli’s equation 
(y—1+2 sin? »)=constant, 


0 (6+ u)/@B=4 (1—du/dth=4 (y+ 1)sec? 
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Moreover, « is constant along the straight Mach line UV, and hence, by (7), 
(1/Hs)y (hay (ha)v] sin 2u= —Nr, (16) 


where r denotes the distance from U to V. Differentiation of (16) gives (since 
dr=h,d2) 


and hence, by (12), (5), and (15), 
H30 08" | — =(1 —N) Ha’, 


and upon integrating we have 


(8’/H3), —(8’/ =(1 —N) | 2’dr, 


ie. by (16), 


L 
8’, + 1—N) | 
G 


where H,° is the value of Hs at G, and r, the distance between G and L (repeated 
use is made of the fact that N is constant on a straight Mach line). The same 
equation holds with f’ and 2’ replaced by 6” and 2”, respectively. But B”;=/'c, 
since the correction does not alter the velocity components on the Mach line BA. 
Moreover, by (11) and (10), 2” has the same value at any point on GL as at the 
point on the axis lying on the same plus Mach line, and on the axis 2” is constant 
and equal to 2’,. Hence, 


L 
B’, | (a’ — 2’,)dr. . (17) 


It is convenient to extend the meaning of r so as to denote by it the distance, 
measured along the minus Mach lines, between the Mach line BA and an arbitrary 
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other plus Mach line in the region ABC (Fig. 1). Then, from Fig. 3, 
V.V,/OV,=U,U,/OU,=cot 2nd (6+ 1). 
since » is constant on a straight Mach line, and hence, 
—rcot 2ud (6+ 4). 
i.e., by (15), ér/dB= —Nrcot2u. 
Moreover, z=constant in the simple-wave region ABC, so 
dB = —2dt= —2du/(1-—2N), 


by (5) and (15), and hence, 


exp | du =f 


where f()= { (sin see } @8) 
r, may be interpreted as the value r takes on AC (Fig. 1) and by x* may be denoted 


the distance from A at which the axis is crossed by the plus Mach line intersecting 
AC at a distance r, from A (Fig. 1). Then 


COS =4xt 
and 
r/f(u)=4x7 sec (m4). 
In particular 


(x, —xe)/[f (4) cos (6+ 2) Je =4 (4) cos (19) 


and (17) may be written 


x 


0 


The correction (14) for the slope of the liner at any point L is seen to consist 
of two parts. The one, 4(2’9- 2’), is due to the difference between the values of 
a on the plus Mach line through L in the uncorrected and the corrected nozzle, 
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respectively. The other is the contribution to the correction at L due to the total 
corrections applied to the liner upsream of Lf. 


The segment of the liner that is placed symmetrically to BC may be treated 
similarly. If / denotes the point with co-ordinates x,=x,, y,=—y,, then the 
correction to this liner segment is 


y*-y"= faa +tan?6,)dx,, (21) 
*B 
and it is found that 
and 
a’, (f (44) C08 |, | (23) 


(use is made of the fact that H. (x,—y)= — Hs (x, y)). 


It may be noted that knowledge of 2™ and 6™ on the axis, in the test rhombus 
alone, is sufficient for the evaluation of the liner corrections (13) and (21). Such 
knowledge is not sufficient for the determination of, for example, the pressure on 
the liner segment BC of the uncorrected nozzle (for the Mach lines of the original 
nozzle design are also the characteristics of the hyperbolic system (8) and (9) and 
hence, the values of 2’ and #’ on the axis in the test rhombus do not determine the 
values of these variables outside the test rhombus). 


5. APPLICATION TO NOZZLE DESIGN 


The most convenient procedure for the evaluation of the liner correction, for 
a nozzle giving a flow of uniform direction before correction, would appear to be 
as follows. 


Consider the upper half of the nozzle (Fig. 1), and let A be the apex of the test 
rhombus, BA and AC the Mach lines through A in the field of flow of the original 
nozzle design, and BC the segment of the liner between them. To any point L on 
this liner segment there corresponds a point Q on the axis lying on the same plus 


tif (6’, —8”, is neglected and x, is determined from x, by the help of-a drawing of the Mach 

line field, Harrop’s®) solution is obtained. If (4’,—8”,) is neglected and L is taken to be the 
point of intersection of the straight line QR with the liner, Ward’s”) solution on Linear Theory 
is obtained. 
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Fig. 4. 


Mach line as L (Fig. 1), and a point G on the Mach line BA lying on the same minus 
Mach line as L. 


The Mach line BA is the last Mach line that has to be computed point by 
point for the original nozzle design, and hence there is a set of points G on it at 
each of which the co-ordinates xg, yg, the Mach angle yu¢, the stream direction 9,, 
the co-ordinates x,, y, of the corresponding point L, and the distance r, between 
G and L have already been computed. For each of these points G find the 
abscissa xg of the corresponding point Q on the axis from (19) and (18). 


If the flow in the test rhombus of the uncorrected nozzle is of uniform direction, 
=constant=, and (by (5)) 2’ =t/+6,. Hence, by (4), 


a’ (dt/dp) (p’e = P's) =- (pq? tan (p™e - (24) 


at any point Q of the axis, where p™ is the pressure in the uncorrected nozzle, and 
p, q, » are the density, velocity magnitude, and Mach angle, respectively in the test 
rhombus of the original nozzle design. Therefore, (8’,—8”,) may be evaluated 
from (20) and (24) when the pressure distribution has been measured on the axis, 
in the test rhombus. In equation (20), x+ denotes distance on the axis, measured 
downstream from A and N is defined by (15), f(«) by (18); Hs must be found by 
numerical differentiation on the Mach line BA if the nozzle has been designed by 
the conventional method of characteristicst. 


T(Hg sin 24)-'=hg=(0x/ 8) sec (6—»), by (6) and (7) where dx/d8 is the ratio of the rates of 
change of x and 8 (=(7/ 180) (¢+P—1000°), where P is Busemann’s pressure number) on the 
Mach line BA. Note that Hs vanishes like /(8—84) near A, if the velocity gradient is 
continuous at A (Appendix of Ref. 5). 


206 


3 
7 
7 
a | 
2 
—-a—--—< -— - — — 
4 
| 
|_| 


CORRECTION OF WIND TUNNEL NOZZLES 


Finally, the correction 5, to the inclination of the liner at L is found from (14), 
and the correction (y* — y”) to the liner ordinate for the abscissa x from (13). The 
opposite correction must be applied to the liner in the lower half of the nozzle. 
The Mach number in the test rhombus of the corrected nozzle will be the same as 
that in the uncorrected nozzle at A. 


6. TEST FOR UNIFORM STREAM DIRECTION 


A check whether the stream direction in the test rhombus of the uncorrected 
nozzle is indeed uniform can be made by measuring the pressure at points in the 
test rhombus placed symmetrically with respect to the axis (Fig. 4). For if 1, 2, 
3, 4 form a Mach quadrangle as shown, then by (5), (10) and (4), 


and in the test rhombus & is equal to the actual stream direction in the uncorrected 


nozzle. (Only differences in & can be measured in this way.) 


If the stream direction is not uniform, the distribution of 6 on the axis may 
be found by pressure measurements, except for an additive constant. For by the 
same argument it also follows that 


and (22) and (23) hold for any four points in the test rhombus such that 1, 2 lie 
on a line parallel to the axis (with 2 downstream of 1) and 3, 4 complete the Mach 
quadrangle (with 1, 4, 2, 3 for the anti-clockwise order of the vertices, if the stream 
comes from the left). 


When &—6, and p™—p™, have been measured on the axis, the corrections 
to the upper liner are found from (20), (14), (13) (and those to the lower liner from 
(23), (22), (21)) as described in Section 5, but instead of (24) the equations 


a’, > 0, = (pq? tan (p™o p™ a) 
and 


which follow from (5) and (4) must be used. 
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SOME DEVELOPMENTS OF EXPANSION METHODS 
FOR SOLVING THE FLUTTER EQUATIONS 


by 


JOHN WILLIAMS, M.Sc. 


(Aerodynamics Division, National Physical Laboratory) 


SUMMARY 


In a theoretical flutter analysis, the real zeros x<0, y>0O, of a stability 
determinant 4 (./x, y) which is a polynomial in y are required, and can readily be 
evaluated once the polynomial expansions are known for assigned values of x. 
When the air loads are estimated on the basis of classical derivative theory, 4 
becomes a bivariate polynomial in / x and y of a special type. Direct expansion of 
the stability determinant becomes too involved when the number of degrees of 
freedom is large and indirect methods of expansion are needed. 


A method of bivariate expansion which uses a framework of lines subject to 
certain restrictions is examined in Part I, and is applied to a quaternary wing-flutter 
problem. The analysis shows that frameworks previously proposed, having all 
intersections outside the flutter quadrant (x<0, y>0O), demand very high 
computational accuracy. A new framework is suggested having all intersections 
inside the flutter quadrant. The improvement in computational accuracy is 
considerable. 


In Part II, methods for the univariate expansion of A with x assigned are 
briefly discussed, and some methods of bivariate expansion are given which provide 
useful alternatives to the framework method. 


1. INTRODUCTION 


In a theoretical flutter analysis, the aircraft structure is usually treated as a semi- 
rigid body with a finite number of degrees of freedom, n say, and the corresponding 
Lagrangian dynamical equations are expressed in non-dimensional form. It then 
follows that the flutter speeds and flutter frequencies of the system are given by the 
real zeros x= —w?<0, y>0, of the n™ order stability determinant 
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where » denotes the non-dimensional frequency parameter and y is a non-dimensional] 
elastic stiffness—air speed parameter’. The elements of the matrices a and e are 
known real constants representing inertias and elastic stiffmesses respectively; 
structural dampings are neglected to avoid further complication of the problem. 
The elements of the matrix .o/(./ x) are in general real transcendental functions of 
x and represent the aerodynamic forces. Sometimes they are tabulated against w, 
but in classical derivative theory they are expressed as quadratics in x as a 
simplifying approximation. 

For convenience, let x=s, X, y=s, Y, where s, and s, are scale parameters which 
may be arbitrarily assigned. The stability determinant may be expressed with 
advantage in the form 

. 
where f(X,Y) and g(X,Y) are real functions of X and Y, and are moreover 
polynomials in Y of degree nm and (m—1) respectively. The real zeros x<0, y>0 of 
A(x, y) are clearly the common zeros X¥<0,Y>0 of the functions f (X,Y) and 
g(X,Y). They can therefore be evaluated readily by examination of the polynomial 
expansions for f(X,Y) and g(X,Y) corresponding to assigned values of X, as 
illustrated in Appendix I. When classical derivative theory is employed, f(X, Y) 
and g(X,Y) become bivariate polynomials of degree nm and (n—1) respectively in 
X and Y jointly, and can therefore be specified completely by the values of 
A(+ at only 4 (m+ 1)(n+2) values of (x, y). 


The problem of obtaining such bivariate expansions and univariate expansions 
(X assigned) is now examined in the case when the number of degrees of freedom 
is so large that direct expansion of A (/ x, y) becomes too involved. 


Notation 


A,,, B,, framework bivariate expansion coefficients corresponding to functions 
f(X,Y), (X,Y). See (5) 


SZ (/x) matrix of non-dimensional aerodynamic coefficients 
a matrix of non-dimensional inertial coefficients 

a,,b, polynomial expansion coefficients for f and g. See (11) and (17) 

@,s,b,, bivariate polynomial expansion coefficients for f and g. See (18) 
C, Lagrangian expansion coefficients. See (14) and (22) 

is sania D, — special expansion coefficients for Frazer framework. See Section 3 
e matrix of non-dimensional elastic stiffness coefficients 
f, g functions related to expanded stability determinant. See (2) 


ft, f- sums of numerically positive and negative terms in expansion for f 
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defined in similar manner to ft, f~, f* 


L, line of framework. See (3) 
n number of degrees of freedom of flutter system 
S,, 8, scale parameters for x and y 


W, operator corresponding to Gregory-Newton formula. See (12), (13), 
(25) and (26) 


X,, Y: co-ordinates of collocation points for expansion by Gregory-Newton 
formula. X,=X,+pv, Y:=Y,+tw 


co-ordinates of framework intersections 
x  —w?*, where w represents a non-dimensional frequency parameter 
y non-dimensional elastic stiffness—air speed parameter 
Z See (10) 
4 - stability determinant. See (1) 
8 operator denoting forward differences 
(X-X,)/v, (Y-Y,)/w 
Il product sign 


= summation sign 


PART I 


A NEW FRAMEWORK OF LINES FOR THE FRAMEWORK METHOD OF 
BIVARIATE EXPANSION 


2. THE FRAMEWORK METHOD OF BIVARIATE EXPANSION 


This method for the bivariate expansion of the stability determinant 
corresponding to classical derivative theory was proposed by Frazer®:*’. Let a 
framework of (m+ 2) lines 

L,=Y —p,X +q,=90, (r=1 tom+2). (3) 


be chosen so that no two lines are parallel and no three are concurrent. Then the 
function f (X,Y) may be expanded in the bivariate partial fractions form 
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or in the bivariate product form 


rs ers 


where k, r, and s range from | to (m+2), and where the 4 (n+ 1) expansion 
coefficients A,, are free for determination. If (X,,,Y,,) denotes the point of 
intersection of the lines L,=0, L,=0, it follows immediately that 


A,; =f (Xr Y,.)/ I L, (Xie (6) 


xr, 


The formule relating to the function g(X, Y) are similar, but only (m+ 1) lines are 
needed and only 4n(n+1) expansion coefficients B,, are involved. For the 
determination of a particular pair of expansion coefficients A,, and B,,, the 
evaluation of two real determinants or of a single complex determinant is required, 
according as X,, is positive or negative. 


The two frameworks proposed by Frazer“: *) and by Price®’, both having all 
intersections outside the flutter quadrant, are discussed in Section 3, and a new 
framework is suggested with all intersections inside the flutter quadrant (i.e. 
X+s<0, Y,,>0). Results obtained by the application of these frameworks to a 
quaternary wing-flutter problem are compared in Section 4. 


One possible method for the determination of the required common zeros 
is by substitution of trial values for (X,Y) in the expansions (4) or (5), but 
speedier methods may be devised in which trial values are assigned only to one 
variable and the remaining variable is eliminated to give a test function. A method 
particularly suited for use with the new framework will be described later. 


3. TYPES OF FRAMEWORK 


If no framework lines cross the flutter quadrant, expansion in the bivariate 
partial fractions form (4) is most convenient. On the other hand, if some lines 
do cross the flutter quadrant, the product form (5) is best as large values occur 
with (4), and trial points can be usefully chosen to lie on these lines. 


The Frazer framework 


The framework proposed by Frazer“) has the properties that all the intersections 
have positive abscissae, and that no lines cross the flutter quadrant. The lines are 
all members of the family 


...-@ 
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where the values of j adopted for the framework are all positive integers and are 
chosen to yield well-separated intersections (X,,, Y,.)*. Frazer usually omits the 
line L»,. Of the framework, in which case the expansions (4) and (5) require slight 
modification and only (n+ 1)? real determinants need be computed. The expansion 
(4) for f (X, Y) is replaced by 


f{(X, L.= = { A,/L,L, } { D,/L, } 


75 


where k, r and s range only from one to (n+ 1). 


The Price framework 

The framework proposed by Price? is composed of lines all crossing the flutter 
quadrant, but having all intersections with positive abscissae. The lines are all 
tangents to the parabola Y?=4X at the parametric points (s? k’, 2 sk) where k takes 
the integral values one to (m+2) and s is a scale parameter which is to be arbitrarily 
assigned*. A useful transformation to new co-ordinates (€,7) may be employed, 
defined by the relations 


A new framework 

It can readily be demonstrated that difficulties regarding accuracy are likely 
to be encountered with both the foregoing frameworks, but that they can be reduced 
by choosing a framework with some of its intersections in the neighbourhood of the 
zeros to be determined. If, in the product expansion (5) for f(X,Y), a trial value 
of (X, Y) is chosen which is close to a zero (X, Y) then f(X, Y) —>0 as 
(X. ¥) —>(X, Y), but the separate terms in (5) will be non-zero and large in general. 
However, these terms could all vanish if one of the framework intersections 
(X,,, ¥+s) happened also to lie at (X,Y). This suggests that if (X,,, Y,.) lay close to 
(X,Y), the terms of the summation would all be small and thus of an order of 
magnitude comparable with the value of f(X, Y). 


The new framework therefore has been chosen to have all its intersections in 
the flutter quadrant, X¥<0 and Y>O. The lines are all members of the familyt 


. . @) 


where the values of j adopted for the framework are all positive integers and are 
selected to yield well-separated intersections (X,,,Y,.). For the treatment of the 
quaternary problem in Section 4, the framework given by the values j=0, 1. 2, 3, 5, 7 
is adopted (see Table I and Fig. 1.). A useful framework of fen lines is given by 
the values j=0, 1, 2, 3, 5, 9, 13, 18, 23, 28. In the product expansion for g(X, Y). 


*The scale parameters s, and s, both have the value unity. 
{This family was obtained from the family used by Frazer by a simple linear transformation. 
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© ZEROS OF STABILITY DETERMINANT 
FOR QUATERNARY PROBLEM 
(s, =I, s,= 30) 
1-0 
Y 
L, 
L, 
=(-§ 


Fig. 1. 
Lines of new framework used for quarternary wing-flutter problem. 


where only (7 +1) of the (n+2) lines used for f(X, Y) are needed, the first line of 
the framework (j=0) should be discarded. For the determination of a particular 
zero of A(/ x.y) in a given problem, s, and s, should be so chosen that the zero 
is likely to occur in the region of the framework intersections. 

In the treatment of the bivariate product expansions obtained for f (X,Y) and 
g(X,Y), a change of variables from (X,Y) to (Z, Y) may usefully be made, by 
means of the relation 


For, when Z is given one of the values adopted for j in the specification of the new 
framework, (X, Y) is effectively assigned to lie on a particular line of the framework, 
L,=0 say; the summation in (5) is then with respect to r only and involves merely 
(n+ 1) terms and expansion coefficients. Let some trial value be assigned to Z. The 
product expansions for f(X,Y) and g(X,Y) may then be arranged as polynomials 
in Y, whose coefficients are linear in the expansion coefficients A,, and B,, 
respectively. The multipliers of A,, and B,, in these combinations depend only on 
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TABLE I 


New framework used for quaternary wing-flutter problem. 


1 0 0 0 
3 2 050 0.05 3 2 02 005 
5 5 125 0.50 4 2 03 0075 
6 7 LIS 105 4 3 04 0.5 
5 04 0 
5 2 O05 0.125 
5 3 06 025 
5 4 07 0375 
=¥ +0.25 j[X +0.1-1)] 
6 2 07 0.175 
6 3 08 035 
6 4 09 0525 
A 


1.1 0.875 


Z and the framework constants, and can be tabulated against Z if the framework 
is standardised. By a reduction scheme involving cross-multiplication between the 
coefficients of the polynomials, a test function T (Z) can be evaluated which vanishes 
when the polynomials have a common real zero, and the corresponding value of 
Y can be simultaneously determined. A computational scheme for the case n=4 
is set out in Appendix I. The relation (10) then gives the corresponding value of X. 


4. NUMERICAL APPLICATIONS 


A quaternary wing-flutter problem was treated in previous calculations with 
the Frazer framework"). In the present paper, the same problem is considered with 
the first six lines of the Price framework (s=1.5), and also with a framework 
consisting of the six lines of Table I, (s.-=1, s,=30). Table II gives the expansion 
coefficients obtained with the three frameworks. 


For a comparison of accuracy between the frameworks a difference factor will 
be introduced. At a chosen (X,Y) value, let f+ and f- denote the magnitudes of 
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the sums of the positive and negative terms respectively in the summation (5) for 
f(X, Y); then the difference factor f is defined as 


f>=[true value of f(X, Y)]/f*. 


If, therefore, a rough estimate of f* is made, it is possible to decide what accuracy 
is needed in the evaluation of f+ and f-, and hence in the expansion coefficients and 
related products, to yield a specific accuracy in f(X,Y). The symbols g*, g~, 2° 
have a similar significance in relation to g(X, Y). 


TABLE HI 


Relative accuracy of expansions given by various frameworks at some trial points 
for quaternary wing-flutter problem 


| f (Frazer) 


x y (Price) f (New) |#(Frazer) | f*(Price) 
-0.1125 | 1.425 |+0.00001 | — |+0.00000152 | 3x10-* | 5x10-™ | 
-045 |12 |-0.00210 | |-0.00210138 | 5x10-> | 5x10-* | 1x10-? 
-0225 | 2.925 |+0.0001 | +.0,00006167 |2x10-7 | 5x10-" | 2x 10-2 
|-0.0069 | ~0,0069839 | 1x10-* | 4x10-* | 1x 10- 
-0.405 | 4.41 | =  |~0,000188297 | 1x10-* | | 610-2 
1.35 42 |—0,003 | - | —0,003796 1x10-* | 1x10-* | 7x 
-0.4455 | 486 |-0.001 - |~0,000960503 | 3x10-? | 3x10-* | 1x10-? 
-1.485 | 4.65 +0,002 | +0.00168 5x10-7 | 4x10-* | 2x10-¢ 

y | g(Frazer)\g(Price)| g (New) g*(Frazer)| g*(Price) | g(New) 
-0.1125 | 1,425 | 0.00031} — 00003150 | 4x10-* | 3x10-7 | 2x10-8 
-0.45 1.2 0.0007 | 0.001 | 0.0006733 1x10-° | 6x10-7 | 2x 10-9 
0,225 2.925 | 0.0025 - 0.0024899 5x10-* | 1x10-* | 1x10-? 
2.7 0.0032 | 0.003 | 0.00325 6x10-* | 1x10-* | 6x 10-4 
-0.405 | 4.41 | 0.002 | 0.002 | 0.00215397 | 1x10-* | 1x10-* | 1x10-! 
-135 | 42 0.007 | 0.006 | 0.00679 4x10-* | 2x10-* | 3x 10-4 
0.4455 | 486 | 0.003 | 0.003 | 0.00274923 | 1x10-* | 9x10-7 | 1x10-? 
-1485 | 4.65 | 0.008 | 0.008 | 0.00792 3x10-* | 2x10-* | 210-4 
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Table III compares values of f=f+—f-, g=g+-—g-, f4, and g* obtained with 
the various frameworks for some typical trial points. It is seen that the accuracy 
in the determination of f (X, Y) and g(X, Y) is much greater with the new framework 
than with the Frazer or Price frameworks. At trial points in the neighbourhood 
of the common zero (x, y)=(— 0.48, 4.86)*, the determination of f (X, Y) and g(X, Y) 
to a particular accuracy requires the retention of at least five more significant figures 
with the Frazer framework and eight more with the Price framework than with the 
new framework. 


5. CONCLUSIONS 


With the new framework, the improvement in computational accuracy is 
considerable. Moreover, the computation required for the treatment of the bivariate 
product expansions can be reduced, since the framework lines cross the flutter 
quadrant and a simple change of variables may usefully be made. The construction 
of single-entry tables, for the general application of this framework to flutter 
problems, is desirable. Tables for the case n=6 are now in preparation and will be 
issued in due course. 


PART II 
OTHER EXPANSION METHODS 


6. UNIVARIATE EXPANSION METHODS 


Univariate expansion in Y is of particular interest in connection with flutter 
analysis at very low frequency parameters, when the approximations of classical 
derivative theory become inadequate. 


When the elements of the aerodynamic matrix are given as tabulated functions 
of w, 4 is polynomial in Y only. If m=<4, A may be expanded directly in Y for any 
assigned X, without much difficulty. For n>4, indirect methods of expansion are 
preferable. With X assigned, write 


r=0 


Then the (n+1) coefficients a, can be calculated by collocation of the function 
f (X, Y) at (n+ 1) distinct points Y=Y,, where t ranges from 0 to n. 


*This was the most important zero from the flutter standpoint. The other common zeros lay 
at (x, y)=(—0.99, 3.57) and (—0.17, 1.38). 
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For example, with equally-spaced collocation points Y,=Y,+tw, the Gregory- 
Newton interpolation formula may be employed. With X assigned, 


where W,,(y) denotes the differencing operator* defined by 


Since »=(Y —Y,)/w, the products in (13) can be easily expressed as polynomials 
in Y. 


With collocation points not necessarily equally spaced (11) may be obtained 
in the Lagrangian form 


(X.Y=2(C, (Y-Y,)} P (14) 
k=0 tk 
tk 


Alternatively, a method of undetermined coefficients may be applied, to provide 
a set of (n+ 1) linear simultaneous equations 


0 


for the coefficients a,. 


Expansion by the Gregory-Newton formula is the fastest of the preceding 
three methods. The formule appropriate to 


r=0 


are similar but only n collocation points Y, are needed. Each of the methods 


*For convenience, the forward difference notation used here is 
(X, Yeu)—f (X, Ys) 


5” F(X, Y:)=8"-'f (X, Yt41)—8"-1f (X, 
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described requires the evaluation of (n+ 1) complex determinants A(X, Y,) for 
the determination of the polynomials f(X, Y) and g(X, Y) at any assigned value of 
X<0. From a computational standpoint, it would be helpful to standardise the set 
of collocation points to be employed, for the products in (13) and (14) could then 
be expressed as polynomials in Y, once for all, and the work required for the solution 
of the equations (16) could be reduced considerably. 


More complicated methods which do not employ collocation can also 
be used, and these are possibly quicker. The stability determinant for an 
assigned x is first reduced to the simpler well-known form |U (x)—yl| where 
U (x)=e-' [—ax—.A(/x)], e~' is the reciprocal of e, and / is the unit diagonal 
matrix. The modified determinant can then be expressed in y directly by methods 
reviewed and discussed by Hotelling’’’ and Wayland". 


7. BIVARIATE EXPANSION METHODS 


When classical derivative theory is employed to give the aerodynamic forces, 
the function f (X, Y) becomes a bivariate polynomial of degree n, jointly in X and Y. 
The coefficients in the Y-polynomial (11) are in fact polynomials of the form 


@,=a,(X)=3 a, X=", . (8) 
s=0 


and in all there are 4 (m+ 1)(m+2) disposable coefficients a,,. The function f(X, Y) 
may therefore be determined by collocation at 4 (n+ 1)(n+2) distinct points (X,, Y,), 
where p ranges from 0 to n but ¢t only from 0 to (n— p). 


Expansion based on univariate interpolation formulae 


With X=X,, the (n+ 1) coefficients a,(X,) in the Y-polynomial for f(X,, Y) 
can be deduced from a Gregory-Newton relation of the type (12) namely: — 


r= 


by comparison of the coefficients of the successive powers of Y in the left and right- 
hand sides. In particular a,(X) — which is independent of X — can be evaluated 
immediately as a,(X,). Next, with X=X,, the relation 


a, ¥)—a, (X,)¥* } y=y, + (20) 


gives the value of a,(X,). Hence, a,(X)=a,,X+a,, is known at both X, and 
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X,, and is therefore completely determinable. More generally, with X =X, at the 
(p+ 1)™ stage the relation 


n 
x a,(X,) — Y./w) f(X,, (21) 
r=0 


gives the values of a,(X,) for r=p to n*. The polynomial a, (X) is now known at 
(p +1) values of X and may therefore be determined completely by collocation in X 
at these pointst. A simple computational scheme for the application of the 
foregoing method is outlined in Appendix II. 


Treatments based on the Lagrangian form (14) or on undetermined coefficients 
can also be devised. At the (p+ 1)" stage the Lagrangian form provides the relation 


=x a,(X,) Y" {Cy (22) 


r=p 
where k and ¢ range from 0 to (n— p), and where 


1 
r= 


Alternatively, by the method of undetermined coefficients, the (n—p+1) linear 
simultaneous equations (t=0 to n—p) are obtained 


n 
r=p 


r 


for the coefficients a, (X,). 
Expansion by the Gregory-Newton formula is the fastest of the preceding three 
methods. 


Expansion based on bivariate interpolation formulae 


With the collocation points X¥,=X,+pv, Y:=Y,+tw equally spaced in the 
X and Y variates separately, an extended form of the Gregory-Newton interpolation 
formula could be employed. It is easily shown? that 


f(X, Y)=W,, (E. n) f (Xo. Yo) (25) 


where £=(X-—X,)/v, n=(Y-Y,)/w, and W,(€,n) represents the differencing 


*Only (n—p+1) collocation points Y=Y+, t=0 to (n—p), are employed at the (p+ 1)™ stage. 


+Alternatively, the value of ap (X) at any other value of X may be determined by interpolation. 
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operator* defined by 


€(€-1) 


Wi m= 1+ + 


n: 


The products involving £ and 7 can readily be expressed as bivariate polynomials in 
X and Y, and the bivariate differences can be built up simply and quickly by 
successive double-entry tabulations, beginning with the values for f (X,, Y,). 


General application 

The formulae appropriate to g(X,Y) are similar to those given for f(X, Y) 
but only 4n(n+1) collocation points are needed. Each of the methods described 
requires the evaluation of 4(m+1)(m+2) complex determinants when all the 
collocation points have negative abscissae. It would be helpful to standardise the 
set of collocation points to be employed, for the products in the Gregory-Newton 
formulae and in the Lagrangian forms could then be expressed as polynomials, once 
for all. A useful set would be 


with the points (X,, Y:) discarded for the determination of g(X,Y). If the values 
of the common zeros to be determined can be roughly estimated beforehand, the 
scale parameters s, and s, should be chosen so that some of the collocation points lie 
in their vicinity, as the computational accuracy will thereby be improved. 


8. CONCLUSIONS 


The bivariate expansion methods given in Part II involve more computational 
operations than the framework method, but the extra work is small compared with 
the amount needed for the evaluation of the determinants required in both cases. 


*This bivariate forward difference notation, 
51° f (Xp, ¥t)—f (Xp, Yo), 
8.1 f (Xp, ¥t)=f (Xp, (Xp, 
™ (Xp, ™ f (X psi, ™ (Xp, Yo), 
= ™1 (Xp, ™-1 (Xp, 


is a natural extension of the difference notation used earlier. 
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Some improvement in accuracy compared with that of the framework method is 
likely, as the collocation points can be distributed evenly over the region in which 
the required zeros are expected to occur, whereas the line intersections of a 
framework are of necessity unevenly spaced. 


6) 
in 
is APPENDIX I 
Cross-multiplication treatment of polynomial expansions 
(n=4) 
When the polynomial expansions 
4) 
n 
“e have a common zero at the assigned X (or Z), the Y-eliminant J (=f, 2g, —f, 2) 
os derived later vanishes, and the corresponding Y-value is given by either of the two 
identical expressions —f,/f,, —g,/g. For the determination of the common zeros, 
it is more accurate to examine the test function T=J/(f,g,) rather than J. 
Alternatively, the vaiues of Y given by the above two expressions could be compared. 
le Checks* 
r=0 ' | 2 3 Operati = 
Sum of Cross-multiplied 
coeffs. sums 
| a, a4 @ @ a, > a, 
h bn ' b DB, = b 
d, d, d, | Di c¢=- & b, >» d, b, = Cr — Cy = b, 
e &.-&¢ c d, d, 
fo fh d,e, — e d, x - eid, 
J f, 8, f, 
*Checks on the numerical working are provided by opposite entries in these two columns. 
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APPENDIX II 


Computational scheme for bivariate expansion method based on 
umvaniate Gregory-Newton interpolation formula 


The scheme to be outlined here involves essentially the formation of simple 
difference tables and some inverse interpolation. The values of f(X,, Y:), where 
p ranges from 0 to n but ¢ only from 0 to (n—p), are assumed known at the outset. 


For the first stage with X¥=X,, the values f(X,, Y;:) are listed against Y,, 
t=0 to n, and the corresponding differences are formed. The coefficients a, (X,) 
for r=0 to n then follow immediately from the relation (19)t. It is convenient next 
to prepare to tabulate each coefficient a, (X) at the values X= X,, p=0 to n, and to 
make difference tables as illustrated below. The column a,(X) can be completed 
immediately, since a, (X)=constant =a, (X,). 


a, (X) | a, (X), a, (X)| | a, (X), a, (X), a, (X) | 


| 


| 


| 


| 


| 
| 
| 


constant 

* 

constant 


The values f(X,, Y:)—a, (X,) Y.", t=0 to (n— 1), can now be listed against Y, 
and the differences evaluated, so that the values of a,(X,) for r=1 to n follow 
immediately from the relation (20). Now the first order differences are constant 
since a, (X) is linear in X, so from the two known values a,(X,) and a, (X,) the 
column a, (X) may be completed by inverse interpolation. 


tIt is desirable that a standard set of collocation points should be employed for which the 
expressions [ 7 (7—1)...(7—m-+1)]/m! have been expressed once for all as polynomials in Y. 
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Similarly, with X = X, at the (p + 1) stage, the values of 


{(Xp¥,) a,(X,)¥,""", t=0 to (n—p), 
may be listed against Y, and the differences evaluated, so that the values of a,(X,) 
for r=p to n follow immediately from the relation (21). But the p™ order differences 
are constant, since a, (X) is a polynomial of order p in X, so from the known values 
a, (X,), @(X,), .... @ (X,), the column a,(X) may be completed by inverse 
interpolation. 


Thus on completion of the (n+ 1)" stage with X = X,, there is to hand a complete 
table of each coefficient a,(X) at the values X=X,, p=0 to n, together with the 
corresponding differences. Hence, the value of each coefficient a,(X) at any 
intermediate X may be found by interpolation, or the corresponding polynomial in 
X may be determined directly, by means of the Gregory-Newton formula. 
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ERRATA 


ERRATUM—Volume II, Part II, August 1950. 


The following correction should be made in Volume II, Part II of the August 
1950 number to “Critical Mach Numbers for Swept-Back Wings,” by S$. Neumark, 
equation (18), page 93: 


For(1=£ 1-sinz) 1 +sinz) 
1+€ 1+sin¢ we € 1-—sin¢g 


i 
\ 
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